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Abstract 

We present results demonstrating the occurrence of changes in the 
collective dynamics of a Hamiltonian system which describes a confined 
microplasma characterized by long-range Coulomb interactions. In its 
lower energy regime, we first detect macroscopically, the transition from 
a "crystalline-like" to a "liquid-like" behavior, which we call the "melting 
transition". We then proceed to study this transition using a microscopic 
chaos indicator called the Smaller Alignment Index (SALI), which utilizes 
two deviation vectors in the tangent dynamics of the flow and is nearly 
constant for ordered (quasi-periodic) orbits, while it decays exponentially 
to zero for chaotic orbits as exp( — (Ai — A2)t), where Ai > A2 > are 
the two largest Lyapunov exponents. During the "melting phase", SALI 
exhibits a peculiar, stair-like decay to zero, reminiscent of "sticky" orbits 
of Hamiltonian systems near the boundaries of resonance islands. This 
alerts us to the importance of the AA = Ai — A2 variations in that regime 
and helps us identify the energy range over which "melting" occurs as a 
multi-stage diffusion process through weakly chaotic layers in the phase 
space of the microplasma. Additional evidence supporting further the 
above findings is given by examining the GALIk indices, which generalize 
SALI (=GALl2) to the case of k > 2 deviation vectors and depend on the 
complete spectrum of Lyapunov exponents of the tangent flow about the 
reference orbit. 



1 Introduction 

It has long been established that microscopic deterministic chaos provides an 
efficient mechanism for the mixing of orbits in the phase-space of dynamical 
systems, leading to the decay of statistical correlations as time evolves. Thus, 
chaotic dynamics can magnify small scale fluctuations and justify the existence 
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of macroscopic variables like entropy and temperature, which are of central 
importance in an analysis based on non-equilibrium statistical mechanics [21] . 

The main purpose of the work presented here is to study the "melting tran- 
sition" in a microplasma model, using standard methods, as well as certain 
recently-developed techniques for chaos detection, such as the Smaller Align- 
ment Index (SALI) (H, l37l, l38l| an d its extension to the so-called Generalized 
Alignment Index (GALI) (H, H^, Q]]. This latter approach is based on geo- 
metrical aspects of the microscopic dynamics and has significant computational 
advantages over more classical indicators based either on local dynamics, such 
as the Lyapunov exponents, or statistical properties, such as the mean temper- 
ature or the Kolmogorov-Sinai entropy. In fact, the use of these novel indices 
often constitutes an improvement of several orders of magnitude in CPU and 
dynamical time for the identification of the chaotic or ordered nature of single 
orbits [HI, [H, [13, 0, 0, & S, ED] . 



The utility of SALI (or GALI) methods in detecting dynamical regime 
changes in few-particle Hamiltonian systems [36[ is due to their sensitivity in 
tracing out the geometrical properties of the tangent dynamics of the flow. 
Thus, they provide accurate information about regime changes when important 
parameters of the system are varied. For example, the existence of "sticky" 
regions and the occurrence of slow diffusion in weakly chaotic domains is de- 
tectable and the distinction between order and strongly chaotic motion is easily 
made by these methods [36]. In the case of fully developed chaotic motion, SALI 
is particularly efficient, since it decays exponentially as exp(— (Ai — \2)t) where 
Ai > A2 > are the two largest Lyapunov exponents [38], while for ordered 
orbits SALI ex const. > 0. 

The system we consider here consists of TV interacting particles described by 
a Hamiltonian function of the form 

H(q,p)=K(p) + V(q) 

where the kinetic energy part K(p) = \ J2i=i Pi 18 quadratic in the general- 
ized momenta p = (pi, . . . ,£>at) and the potential energy V(q) is a function of 
its generalized position coordinates q = (qi, . . . , qw)- In the case of confined 
systems, one takes V as the sum of two terms: V(q) = V tr (q) + Vi n (q), where 
Vtr(q) represents the potential of the trap (being a smooth positive function), 
while Vi n (q) accounts for the interactions amongst the N particles. The fact 
that TV is a finite number classifies the system as "small", in contrast to a "large" 
thermodynamic system (where one lets N and the volume V tend to infinity in 
such a way that N/V is a finite real constant) [24l |23|. 

Small, finite Hamiltonian systems in contrast to large, infinite, systems do 
not exhibit phase transitions in the standard sense. The characteristic disconti- 
nuities or singularities in thermodynamic functions or their derivatives, which is 
the signature of a phase transition in large systems would now appear as steep 
but continuous changes in their thermodynamic functions or their derivatives 

& 

Phase transitions in small systems, have also been associated, quite recently, 
with certain topology changes in configuration space 13, These devel- 



opments are stemming from earlier work on Hamiltonian systems exhibiting 
chaotic instabilities associated with singular behaviour of their configuration- 
space curvature fluctuations at their phase transition point 13, H, 0|. Fur- 



thermore, investigations of the temperature dependence of the largest Lyapunov 
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exponent and other observables related to the "topological hypothesis" and the 
issue of phase transitions in many-degrees of freedom systems is thoroughly 
presented in [l6[ and references therein. 

In this paper, our model Hamiltonian system describes a microplasma char- 
acterized by long range (non shielded) Coulomb interactions described by the 
potential Vi n (q). The microplasma is confined in a Penning trap given by the 
potential V tr {q) and our aim is to study the motion of its ions as the energy 
increases. The system evolves, in general, in 3 spatial dimensions and contains a 
relatively small number of ions for which the tracking of individual trajectories 
is numerically as well as experimentally feasible [H, 26L [ToL fl3L l33l I2H] . 

It is important to note that Lyapunov exponents have already been used 
to investigate spatially extended pla smas [5|, |4|, |42j as well as 1-dimensional 
wave-particle plasma models [13, Moreover, as Gaspard [Hf has demon- 
strated, for a realistic Hamiltonian model in 3 spatial dimensions containing 
a relatively small number of ions, the long range nature of the Coulomb in- 
teraction Vi n (q) makes the maximum Lyapunov exponents behave differently 
in microplasmas than in many-particle systems with short-range interactions, 
such as a hard-ball fluid or systems with shielded ionic interactions (Yukawa 
potentials). A detailed analysis of the differences between shielded (Yukawa- 
like) and unshielded Coulomb potentials can be found in the book by Baus and 
Tejero [6|]. Indeed, in microplasmas composed of more than a few dozen of ions, 
the behavior of the maximum Lyapunov exponent, especially in the "gas phase", 
was explained following purely statistical mechanics arguments (22| . 

The plan of the paper is the following: The Hamiltonian representation of 
our system and a summary of previous work are presented in Sec. 2. Some 
background material related to the SALI method and the spectrum of Lya- 
punov exponents are given in Sec. 3. The calculation of Lyapunov exponents 
spectra for our system is presented in Sec. 4, while the detection of the "melting 
transition" as a passage from weak to strong chaos is demonstrated in Sec. 5. 
Finally, our conclusions are presented in Sec. 6. 



2 Description of the model and summary of pre- 
vious work 

Let us consider a microplasma of N ions of equal mass m = 1 and electric charge 
q in a Penning trap with electrostatic potential 



2 Z 2 — x 2 — v 2 

<S>(x,y,z)=V „ 2 ^ 2 / (1) 



and constant magnetic field along the z direction with a vector potential of the 
form 

A(x,y,z) = ^(-By,Bx,0). (2) 
Then, the Hamiltonian of the full system is given by 



n = f2{^(Pi-<lA(r i )f)+ q $(rA+ £ 

i=l I J l<i<jt 



<j<N 47re °^' 



Q 2 

(3) 
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where n is the position of the ith ion, is the Euclidean distance between the 
ith and jth ions and eo is the vacuum permittivity. In the Penning trap, the 
ions are subjected to a harmonic confinement in the z direction with frequency 



^ (4) 



while in the perpendicular direction (due to the cyclotron motion) they rotate 
with frequency uo c = qB/m. Thus, in a frame rotating around the z axis at the 
Larmor frequency uol = cj c /2, the ions feel a harmonic confinement of frequency 

u x = u y = ^4 — if in the direction perpendicular to the magnetic field. In 

the rescaled time r = co c t, position R = r/a and energy H = m ^i a 2 with 
i 

a = ^ 47re ^ ma;2 ^ 3 , the Hamiltonian (|3]) describing the motion takes the form 

N -i N -i 2 2 i 

i=l i=l i<j J 

where E is the total energy of the system, = (Xi,Yi, Zi) and = {Pxi •> ^Yi > Pzi ) 
are the positions and canonically conjugate momenta respectively, Rij is the Eu- 
clidean distance between different ions i,j given by 



Rij = ^{Xi - XjY + (Yi - Yjf + (Zi - Zjf (6) 

and 7 = uo z /uo c . 

The ions are trapped in bounded motion under the condition that 

< |7l < ^- (7) 

The trap is called prolate if < |7| < isotropic if |7| = and oblate if ^= < 
|7| < So, the motion is quasi 1-dimensional in the limit 7^0 and quasi 
2-dimensional in the limit 7 — > 1/ y/2. The Z direction is a symmetry axis and 
hence the Z component of the angular momentum Lz = Y^iLi XiPyi — YiPxi 
is conserved, being thus, a second integral of the motion. We suppose from now 
on that the angular momentum is equal to zero (i.e. Lz = 0) and that the 
motion is studied in the Larmor rotating frame. 

In [lH, it was shown that the motion of the ions governed by Hamiltonian 
([5]) is generally very sensitive to initial conditions and has at least one posi- 
tive maximum Lyapunov exponent, while a first study of its dependence on the 
energy, number of ions and trap geometry was also presented. At low kinetic 
energies, where the microplasma forms an ion crystal, chaos is considerably 
reduced, since the motion is quasi-harmonic around stable equilibrium config- 
urations. On the other hand, at high temperatures (or kinetic energies), the 
maximum Lyapunov exponent decreases, as Coulomb interactions become neg- 
ligible and the microplasma forms a thermal cloud of nearly independent ions 
moving throughout the full extent of the harmonic potential of the trap. Fi- 
nally, for intermediate values of the energy, there is a regime of chaotic behavior 
which becomes wider as the number of ions increases. It is in this regime that 
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the maximum Lyapunov exponent attains a peak, whose value increases as a 
function of the number of ions. 

A crucial aspect related to our work is that of the geometry of the trap. 
As has been demonstrated in [22] for prolate traps the maximum Lyapunov 
exponent (and therefore the K-S entropy of the system) show a distinctively 
smoother increase in comparison to the oblate traps in the energy values interval 
from E min (its minimum) to Eq (its maximum). The transition of interest for 
our analysis is taking place within a small interval of energies at the beginning, 
i.e. just after E min . Therefore, the slower the increase of the function E = E{H) 
the higher the resolution required and the more accurate (and time-consuming) 
the numerical computation needed. Of course, if one argues in analogy to the 
standard phase transition theory, the dimensionality of the problem makes it 
easier to detect transitions in two or three spatial dimensions. But systems with 
long range interactions may still exhibit a transition even in one dimension. 
This justifies the choice of small 7 adopted in the present paper as being both 
physically challenging and numerically tractable. 

2.1 Crossover through different dynamical regimes 

At zero temperature, the system freezes at a crystalline state (see (22] and in 
particular refs. [20] - [25] therein) reminiscent of Wigner crystals [44|]. In that 
case, the ion crystal is composed of several concentric rings for the case of an 
oblate (quasi 2-dimensional) trap or a single line of ions in the case of a prolate 
(quasi 1-dimensional) trap. As temperature rises slightly above zero, the ions 
execute quasi-periodic motion about stable periodic orbits moving around their 
equilibrium position. This is the regime where stable normal modes of vibration 
play an important role, as expected. 

At slightly higher temperatures, the system "bifurcates" from the regime 
of quasi-harmonic motion to one where it executes quasiperiodic motion in 
the form of a collective soft-mode. We use the term "soft-mode" here since 
the dynamics is similar to what would be expected under a collective force on 
the particles, which diminishes as the energy increases. The frequencies of the 
normal vibrational modes become smaller and smaller and tend to zero at the 
transition. 

As the energy continues to increase, these soft modes overcome the energy 
barriers and set out to explore broader domains of phase space. At even higher 
temperatures, their motion becomes erratic and the ion crystal "melts", entering 
a regime of strongly chaotic motion in phase space. 

As the temperature increases further, the ions form a thermal cloud in which 
the mean Coulomb potential energy starts to become negligible with respect to 
the mean kinetic energy and mean harmonic potential energy. It has been shown 
that the maximum Lyapunov exponent Ai reaches its peak value at energies well 
above the "melting phase", while this transition is not manifested in the behavior 
of the Lyapunov exponents as a function of the number of particles TV or the 
energy E. 

Finally, for energies beyond the peak value of the maximum Lyapunov ex- 
ponent, the motion is dominated by the kinetic energy part and the system 
becomes amenable to a statistical mechanical treatment with considerable ac- 
curacy, in spite of the small number of particles present. Indeed, at such high 
temperatures T, the spatial disorder of the microplasma can be described in 
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terms of its thermodynamic entropy and its maximum Lyapunov exponent Ai 
expressed by the theoretical estimate [Hf 



N 2 \ 2 (lnT)2 

\ - N 3 for T -> oc (8) 



Rii / Ti 



which may be used to characterize the crossover to the regime of thermal cloud 
motion, albeit at a high computational cost for the determination of Ai. 
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Figure 1: Plot of typical orbits for N = 5 and 7 = 0.07 (prolate quasi 1- 
dimensional geometry): (a) At E = 2 the motion is dominated by its stable 
normal modes and "crystalline" behavior is observed, (b) At E = 2.35 a tran- 
sition to chaotic motion occurs, where "melting" is expected to take place, (c) 
At E = 9 the onset of "thermal cloud" behavior is evident. 



Furthermore, in [22], the maximum value of Ai is numerically observed to 
obey a power law with respect to the energy as well as the number of ions. 
Thus, one reaches the conclusion that, for a prolate trap geometry (quasi 1- 
dimensional) , the "melting transition" and its associated microscopic dynamical 
behavior should be detectable at low energies and few number of ions (e.g. for 
N = 5 and 7 = 0.07). 

The different characteristic kinds of motion in such a setting are depicted 
in Fig. [TJ Panel (a) presents a typical behavior of the ions at energy E = 2, 
where the dynamics is dominated by stable normal modes that correspond to 
"crystalline" behavior. At the slightly increased energy E = 2.35 of panel (b) 
we observe the occurrence of a transition from stable collective (quasi-periodic) 
motion to chaotic behavior. Finally, in panel (c), the onset of a fully "thermal" 
cloud behavior is apparent at the energy E = 9. So, for the rest of the paper 
we consider the case of N = 5 ions in the prolate trap geometry with 7 = 0.07. 
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Figure 2: Plot of the dependence of the temperature T (as in eq. (|9]), dimen- 
sionless scale) on the energy E for the microplasma system ([5j) of N = 5 ions 
in a prolate Penning trap (7 = 0.07). The error bars represent one standard 
deviation from the time averaged temperature T. Clearly, the crossover from 
ordered to weakly chaotic motion around E G (2, 2.5) does not occur in a man- 
ner justifying the use of statistical mechanics considerations. In the inset, the 
ratio of the standard deviation of the time averaged T to the time averaged 
temperature T shows that they are comparable. 

The temperature, in dimensionless units, 

rcx Vr 2 ms = (l|P|l 2 > w ^El|P# (9) 

i=l 

(where ||P;|| denotes the Euclidean norm of the zth ions' velocity P^) averaged 
over the time t which is a proportional quantity of the mean kinetic energy of 
the ions is plotted in Fig. as a function of the energy E for the microplasma 
([5]) with TV = 5 ions in the prolate Penning trap with 7 = 0.07. Evidently, the 
melting of the crystal here is not associated with a sharp increase of T at some 
critical energy, as would be expected from a first order phase transition. Note 
also in Fig. [2] that the fluctuations of T (measured by its standard deviation) 
cannot provide a clear indication of the crossover from a regime of stable oscil- 
lations, through "melting", to strongly developed chaotic motion. So, the use 
of statistical mechanics considerations must focus on the importance of large 
fluctuations. Indeed, the inset of Fig. [2] shows that the fluctuations of the 
time averaged temperature T, in the crossover regime, are comparable with the 
values of the time averaged temperature itself. 

One can, therefore, say that there is no "macroscopic" methodology for de- 
tecting dynamical regime changes in system ([5j) that might be useful to theo- 
rists or experimentalists. It, therefore, becomes especially important to adopt a 
different approach and attempt to study dynamical regime changes in our mi- 
croplasma system, by performing a detailed study of its microscopic dynamics. 

To this end, we shall focus hereafter on the lower energy range of Hamilto- 
nian ([5]), where the transition from "crystalline-like" to "liquid-like" collective 
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behavior takes place. In the next section, we review the methods used for the 
detection of these dynamical regimes, i.e. the Smaller Alignment Index (SALI) 



35|, |37|, |38] and the spectrum of Lyapunov exponents [7, 8] and compare the 



information provided by these approaches. 

As we shall see, in the crossover regime, the SALI exhibits a rather intricate 
"stair-like" power law decay to zero, with different Lyapunov exponents, rem- 
iniscent of the "stickiness" behavior typically observed in the neighborhood of 
resonance islands of Hamiltonian systems [36j]. In our case, this type of weakly 
chaotic dynamics will turn out to characterize the energy range over which the 
passage from ordered to irregular dynamics occurs, through what we call the 
"melting transition". 



3 Methods for distinguishing between order and 
chaos 

3.1 The Smaller Alignment Index (SALI) 

The SALI method was initially introduced in [35| and has been applied success- 
fully to distinguish between ordered and chaotic dynamics in maps of various 
dimensions (Til [lift, i n Hamiltonian systems OS [lift, as well as in problems of ce- 
lestial mechanics (iol IZlft . galactic dynamics [Hj], field theory [19] and nonlinear 
1-dimensional lattices^, 0, Elft. 

Following (lift , one considers the 27V-dimensional phase space of an arbitrary 
autonomous Hamiltonian system 



H(q,p) = H( qi (t), . . . , q N (t), Pl (t), . . . ,p N {t)) = E 



(10) 



where q is the canonical position vector, p is the corresponding canonical con- 
jugate momentum vector and E is the total constant energy. The time evo- 
lution of an orbit of Hamiltonian ([TO]) associated with an initial condition 
x(to) = (qi(to), . . . ,qN(to),pi(to), . . . ,pN(to)) at time t = to is defined as the 
solution of Hamilton's equations of motion 



dqi(t) dH d Pi (t) 



dt 



dpi(ty dt 



dH 



i= 1, 



(11) 



and is called the orbit x{t) of eqs. ([TT]) passing through x(to). 

To define the Smaller Alignment Index (SALI), one uses the variational 
equations, which represent the linearization of Hamilton's equations of motion 
([TT]) about a reference orbit x{t) of the system and are defined as 



dvj(t) 
dt 



= J(x(t)).Vi(t), Vi = l,...,27V 



(12) 



where J{x{t)) is the Jacobian of the right-hand side of eqs. ([TT]) calculated 
about the orbit x(t). The vectors ^(t),Vi = 1, . . . ,27V are known as deviation 
vectors of the x(t). If we follow two of them, say Vk{t) and vi(i), we can define 
the Smaller Alignment Index (SALI) as 



SALI(t) ee min 



in I 



v k (t) v t (t) 



\m)\\ mm 



v k (t) Vl(t) 



urn wmw 



(13) 
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where || • || denotes the usual Euclidean norm defined in R 2N . 

If x(t) is chaotic then lim^oo SALI(t) = min{0, 2} = since both deviation 
vectors tend to align with the direction of the maximum Lyapunov exponent as 
t increases 3^, E3l - Furthermore, it has been shown that the time evolution of 



the SALI for chaotic orbits tends to zero exponentially at a rate related to the 
difference of the two largest Lyapunov exponents Ai and A2 as (38| 

SALI(t) oce- (Al - A2)£ . (14) 

On the other hand, if the orbit is ordered, SALI exhibits small oscillations 
about a constant a G (0, a/2]. This is so, because both deviation vectors tend 
to become tangential to the torus on which the orbit is evolving while remain- 
ing mutually linearly independent in time due the existence of local or global 
integrals of motion [371 ]. 

It follows that this sharply different behavior of the SALI for ordered and 
chaotic orbits makes it a reliable and computationally fast tool for distinguishing 
between order and chaos in Hamiltonian systems of any number of degrees of 
freedom. The choice of the initial deviation vectors is arbitrary and in general 
does not affect the method apart from some very special cases demonstrated 



theoretically and verified numerically in a number of previous works [35, 37. 



3.2 The spectrum of Lyapunov exponents 

One of the most standard and well-established means for extracting information 
about the nature of a given orbit of a dynamical system is to calculate its 
maximum Lyapunov Characteristic Exponent (LCE) Ai. If Ai > the orbit is 
characterized as chaotic. The theory of Lyapunov exponents was first applied 
to characterize chaotic orbits by Oseledec [30], while the connection between 
Lyapunov exponents and exponential divergence of nearby orbits was given in 
0, [H]. Benettin et al. @, Q studied the problem of the computation of all 
LCEs theoretically and proposed in Q an algorithm for their efficient numerical 
evaluation. In particular, Ai is computed as the limit for t — > 00 of the quantity 

t (0) || t^oo 

where ^1(0), v*i(t) are deviation vectors from the orbit we want to charac- 
terize, at times t = and t > respectively. Similarly, all other LCEs, 
A2,A3, . . . , A27V are computed as limits for t — > 00 of analogous quantities, 
L2(t),Ls(t), . . . , L2iv(£). In the present paper, we shall compute the values of 
all Lyapunov exponents (called the Lyapunov spectrum) using the algorithm 
proposed by Benettin et al. @, Hj]. 

Let us also recall that the Lyapunov spectrum of an N degree of freedom 
Hamiltonian system, i.e. {K}i=i, exhibits a basic symmetry with N — 1 of its 
members having the opposite sign of the other N — 1 and two exponents, i.e. 
A at, Ajv+i, being equal to 0. So, the discussion from here on concerns only the 
positive half of the Lyapunov spectrum, {A^}^ 1 where A^ > 0. 

The calculation of the Lyapunov spectrum is computationally demanding, 
as its convergence often requires the integration of trajectories over very long 
time intervals. This convergence depends on the inverse of the largest Lyapunov 
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exponent or Lyapunov time. So, the closer Ai is to zero the longer the integration 
is needed to obtain reliable estimates for the full spectrum. 



On the other hand, of particular importance to statistical mechanics is the 
connection of the sum of all positive Lyapunov exponents (an index of chaos 
based on microscopic quantities) to the Kolmogorov-Sinai entropy Hks (a mea- 
sure of disorder of a macroscopic nature), given by Pesin's celebrated theorem 



H KS = Yl A * ( 16 ) 

Ai>0 



in the case of zero escape rates. This connection has opened new directions in 
the applications of thermodynamics, u sing the underlying microscopic dynamics 



of chaotic orbits (for a review see [18l l2ll|). 



4 Calculation of Lyapunov spectra 



For the purposes of simplifying our analysis, let us consider the case of few ions, 
e.g. N = 5 and 7 = 0.07 (i.e. a small system in a prolate trap). As mentioned 
above, in this setting, it is easier to detect dynamical regime changes, such as 
seen in Fig. 02 where we have plotted the Kolmogorov-Sinai entropy Hks (see 
eq. ([T6|) ) as a function of the energy E of the microplasma Hamiltonian ([5]). 
Both Hks and Lyapunov exponents grow less steeply than in the oblate case 
(quasi 2-dimensional trap) for the region of small energies, as already pointed 
out in 

This slow increase of Hks and A^'s exhibits an inflection point at low energies 
as is evident in the inset of Fig. 03 We remark that this inflection point occurs 
just above the energy threshold that permits ions to move around their fixed 
points, while beyond this energy we observe the transition from weak to strong 
chaotic behavior. Note that Hks as well as the Lyapunov exponents exhibit a 
maximum around E w 6.2. 



10 




Figure 3: The Kolmogorov-Sinai entropy Hks, for 7 = 0.07 and N = 5, com- 
puted by the eq. (fl6|) (dashed highest line) and the Lyapunov exponents (black 
solid lines) are plotted from the maximum down to the lowest one (15th) as a 
function of the energy. Note that, as shown in the inset, their maximum occurs 
at energies close to E w 6.2, which is much higher than the regime where the 
"melting transition" takes place (see text). 



As a first attempt to derive useful information from the spectrum of Lya- 
punov exponents (see Fig. [3]), in the energy range where the "melting transition" 
is expected to occur, we calculate the logarithms of the ratios of successive pair- 
wise differences of Lyapunov exponents 



A; - A., 



+1 



Ai+i — Aj+2 



l,...,N-2, Xi = Xi(E) 



(17) 



as a function of the energy (see Fig. [4j). The calculation of all A^'s was carried 
out for each trajectory up to a final integration time tf = 10 6 , ensuring that a 
relative convergence to 4 or 5 significant decimal digits has been achieved. 



Note in Fig. [H(a) that the ratios (fT7|) fluctuate wildly initially and only for 
E > 2.3 do they start to settle down to small oscillations about their ultimate 
values, which, at least for the first 2 or 3 ratios appear to be quite distinct. 
This suggests that it is within an interval Eo < E < 2.3 where Eq — 1-89 that 
we should expect to find the "melting transition", where the positive Lyapunov 
exponents are very small, far from the values they attain in the regime of strong 
chaos (see inset of Fig. [3] for E > 2.3). 
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Figure 4: (a) Plot of the logarithms of the ratios Vi(E), i = 1,...,9 (see 
eq. (|T7|) ). for 7 = 0.07 and N = 5, as a function of the energy E. Note the 
characteristic gaps in the spectrum appearing for large values of E. (b): Same 
as in panel (a) for i = 1, ... ,4 and small values of E, with the characteristic 
gaps in the spectrum becoming erratic in the regime of weak chaos, settling 
down to more clearly defined values around E > 2.3. 



Undoubtedly, although the information concerning the existence of two dis- 
tinct dynamical regimes at small energies is somehow contained in the plots of 
Fig. HI expressions (|T7|) still fail to provide a precise picture of this change. 
This is perhaps due to inaccuracies of the inherent averaging process used in 
the calculation of the Vi{E) at low energies, where the Lyapunov exponents are 
very small. We, therefore, proceed to study more carefully the underlying local 
microscopic dynamics by employing the SALI method we introduced previously. 



5 Weak chaos detection 

In this section, we apply the procedure outlined in Sec. 3.1 to analyze micro- 
scopically trajectories of the few-particle (i.e. N = 5 ions) microplasma system 
described by the Hamiltonian ([5]) in the prolate trap with 7 = 0.07. We are 
thus able to identify at low energies a regime of weak chaos, through which the 
system passes from a "crystalline-like" type of ordered motion to the strongly 
chaotic behavior of collective soft modes. This is achieved by detecting quali- 
tative changes in the SALI behavior for trajectories within the energy interval 
where the transition occurs. 

Fig. [5] shows the behavior of the SALI as a function of time t at selected 
number of energies. These plots represent three distinct classes of long time 
behaviors for typical trajectories within the interval 1.9 < E < 3.5. We have 
also computed, for the same sample of energies, the corresponding Lyapunov 
spectra for sufficiently long integration times (typically tf = 10 6 ) to make sure 
that the Lyapunov exponents have converged to their limiting values up to a 
desired accuracy (typically to 4 or 5 significant decimal digits). 

Fig. [5](a) represents the SALI evolution for a representative trajectory of 
system ([5]) at E = 1.9. We see that it fluctuates around non-zero positive 
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values up to an integration time of order tf = 10 6 , indicating that the motion 
(at least up to that time) is ordered and quasi-periodic. By contrast, in Fig. 
03(b), we follow an orbit with energy E = 3.5 and find that SALI decays to zero 
exponentially fast, driven by the first gap in the Lyapunov spectrum of Fig. [H 
according to the formula e~^ Xl ~ X2 ^ [38j]. In this case, the first two Lyapunov 
exponents of the trajectory, Ai ~ 0.03112 and A2 ~ 0.01746, are large enough 
and distinct, implying strongly chaotic motion. 

Finally, let us turn our attention to the intermediate and more interesting 
case where SALI exhibits a stair-like decay to zero as a function of time, at 
energies 2 < E < 2.3 (see Fig. G3(c)). Choosing, for example, a trajectory with 
E = 2.033 we observe that SALI differs significantly compared to what is shown 
in Fig. [5](a) or 03(b). At first, one might think that, as in Fig. [5](a), the orbit 
is quasi-periodic, since SALI "sticks" to non-zero values for a fairly long time 
interval (t < 4.4 x 10 5 ). Then, its behavior changes qualitatively, showing a near 
exponential decay to zero oc e~( A i- A 2)^ This suggests that, even though the 
Lyapunov exponents (and their differences) are very small, eq. ([T4|) still holds, 
as in the strongly chaotic domains of phase space [381 ]. 
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Figure 5: Plot of the evolution of SALI as a function of time t for TV = 5, 
7 = 0.07 and 5 typical energies: (a) E = 1.9 where the motion is ordered, (b) 
E = 3.5 (strong chaotic behavior) where we also plot the theoretical prediction 
SALI(t) ex e -( A i- A 2)* for comparison, (c) E = 2.033, E = 2.05 and E = 
2.283, plotting also the corresponding theoretical predictions. Note that all 
three vertical axes are logarithmic. 
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If we slightly increase the energy, however, to E = 2.05, SALI displays an 
intermediate and rather intricate behavior, following a stair-like decay to zero, 
shown by the middle curve in Fig. HJc). Here, the exponential law (|T4|) does 
not explain the SALI decay to zero. It appears that the motion lies within a 
weakly chaotic domain and "sticks" temporarily to islands of regular motion, 
executing a multi-stage diffusion process [36j]. Finally, by increasing the energy 
to E = 2.283, we observe that SALI decays to zero exponentially fast, following 
closely the theoretical prediction e^ 1- * 2 ^, with Ai « 0.00166 and A 2 « 0.0007. 
These results imply that the corresponding trajectory is fully chaotic and the 
system has "melted", passing to a regime with strongly mixing properties. 

We now examine more closely this step-wise decay of the SALI, observed 
in Fig. [5](c), over a range of energies spanning the melting transition. To this 
end, we start with our microplasma in the form of a Wigner crystal, at Eq and 
proceed to the onset of thermal cloud formation, increasing the energy by steps 
of AE = 0.05 up to E max 6.2. 




Figure 6: Plot of the quantity Q = for 7 = 0.07 and TV = 5, as a function 
of the energy E G where AA is the difference of the two largest Lya- 

punov exponents and K is the slope calculated by linear regression performed 
on the log(SALI(t)) vs. t curves. We also plot the line Q = 1 to indicate the 
energy range where AA = Ai — A2 ~ K and eq. (fl4]) are satisfied. 



One way to quantify the departure of the SALI from the exponential decay 
law ([T4|) (see Fig. G3(c)) is to calculate for each energy in the interval [Eq, E max ] 
the quantity 

Q = f (is) 

defined as the ratio of the difference between the two largest Lyapunov exponents 
AA = Ai — A2 and K, the linear regression estimate of this difference, obtained 
from semi-logarithmic plots of SALI, like those depicted in Fig. \Mb) and (c). 
The result is shown in Fig. E3 As expected, for the strongly chaotic cases where 
K w AA, we obtain values of Q near 1 with clearly negative log(SALI) slopes. 

The interesting part of Fig. O however, is observed in the energy interval 
2.0 < E < 2.5 where K becomes one order of magnitude smaller than AA. 
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It is in this regime that we observe the "melting transition" for our system. 
The divergence of SALI from its theoretical estimate of eq. ([T4|) signifies the 
presence of other collective modes of motion. This motion is organized around 
islands of stability, suggesting a certain type of "stickiness" around different 
tori at successive time intervals. During these time intervals SALI decays in a 
step-like manner by different power laws. 



In fact, we can provide more evidence to support the above interpretation of 
the dynamics in the "melting phase", using a recent generalization of the SALI, 
called the Generalized Alignment Indices (GALI) [39j]. These indices represent 
the volumes of a parallelepiped, called GALI&, formed by k > 2 unit deviation 
vectors, which (a) for ordered motion, oscillate about positive constants as long 
as k < d (d being the dimension of the torus) and decay following power laws 
for k > d, while (b) for chaotic orbits, they all vanish exponentially with expo- 
nents that become more negative as k increases and depend on more Lyapunov 
exponents (39I. l36j|. In practise, since the computation of the GALI requires the 
calculation of a fairly big number of large determinants at every time step, we 
have adopted an alternative way of evaluating it which is significantly faster in 
CPU time, the so called Linear Dependence Index (LDI) [lj]. 



Thus, in Fig. we use this approach to compute the GALI/e, k = 2, . . . , 10 
and find that they demonstrate the existence of an (at least 6-dimensional) torus 
for E = 1.9 [36j], while they show a step-wise decrease for E = 2, which continues 
to hold up to E < 2.283, as predicted by the SALI (= GALI2) calculations 
(compare Fig. [U(a), (c) with Fig. [T^a) and (b)). Of course, as the energy 
increases, the (weakly) chaotic nature of the orbits becomes evident at earlier 
times (e.g. for E = 2 we have t« 3 x 10 6 and for E = 2.283 we get t« 6 x 10 4 
for the LDI (or GALI) to become ~ 10 -11 ) as we see in Fig. [3(b) and (c). 
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Figure 7: Plot of the LDI/e (=GALI/c), k = 2, . . . , 10 as a function of time t for 
the microplasma system ([5]) of N = 5 ions in a prolate trap (7 = 0.07) for 3 
typical energies, (a) LDI2-LDI10 for the energy E = 1.9. (b) Same as in panel 
(a) for the slightly bigger energy E = 2. (c) Same as in panel (a) for the bigger 
energy E = 2.283. Note that all axes are logarithmic. 

Thus, the above results suggest that the behavior of the SALI during dynam- 
ical regime changes is especially useful, as it alerts us to examine quantities like 
Q (see eq. (fl8|) ) . related to variations in the differences of Lyapunov exponents. 
Such quantities indeed furnish important information, not readily available by 
other more standard microscopic chaotic indicators or by related macroscopic 
entropic quantities. 

6 Conclusions 

In this paper, we have reported results demonstrating the occurrence of dynam- 
ical regime changes in a Hamiltonian system describing a microplasma confined 
in a prolate quasi 1-dimensional configuration and characterized by long-range 
Coulomb interactions. More specifically, in the lower energy regime, we macro- 



16 



scopically detected the transition from "crystalline-like" to "liquid-like" behav- 
ior, through what we call the "melting phase". As expected, well beyond this 
phase, the microplasma exhibits strong chaotic behavior that may be described 
by the macroscopic variables of statistical mechanics. The question, therefore, 
is how can one determine the dynamical nature and identify the energy range 
of this "melting transition" from ordered to strongly chaotic behavior. 

To this end, we first showed that the "melting" of our quasi 1-dimensional 
crystal (composed of few ions confined in a prolate trap) is not associated with 
a sharp increase of the temperature at some critical energy, as might have been 
expected. Furthermore, the positive Lyapunov exponents (and the Kolmogorov- 
Sinai entropy expressed by their sum) attain their highest values at energies 
much higher than the regime where the microplasma "melts". Thus, it appears 
that there is no clear "macroscopic" methodology for identifying and studying 
this "melting" process in detail. 

Next, we argued that even though information about dynamical regime 
changes must be contained in the Lyapunov exponents, Lyapunov spectra by 
themselves still fail to provide, a more precise picture of this change. Indeed, 
while the ratios of the larger of them show strong fluctuations at low energies, 
these require further analysis before they can reveal, with any precision, the 
energies over which the "melting transition" occurs. 

For these reasons, we found it useful to employ a more accurate tool, called 
the SALI method, to study in more detail the local microscopic dynamics of the 
microplasma system. In this way, we discovered an energy range of weak chaos, 
where the positive Lyapunov exponents are very small and SALI exhibits a stair- 
like decay to zero with varying decay rates. As SALI(t) ex e~( Al_A2 ^ in chaotic 
domains, this inspired us to look more closely at the statistical fluctuations 
of the differences of the two largest Lyapunov exponents AA = Ai — A2 in 
that regime. We, thus, observed that these differences exhibit their largest 
fluctuations over a definite energy interval, where "melting" occurs in a way 
reminiscent of "sticky" orbits, executing a multi-stage diffusion process near the 
boundaries of resonance islands in the phase space of Hamiltonian systems. The 
above results were also supported by the use of an extended set of indices which 
generalize SALI to the case of more than 2 deviation vectors. 

Finally, we remark that it is also the rapid convergence of these indices, 
which turns them into efficient diagnostic tools that may be used to replace 
demanding molecular simulations in identifying weakly chaotic regimes in multi- 
particle systems. As we intend to show in a future publication, these indices 
can provide a viable alternative to the notoriously time consuming simulations 
required to study the presence of weak chaos and slow diffusive effects in the 
dynamics of met ast able states. 
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